% Scripts that read res files, assemble result matrix ifft the date back to
% time series
%simulation for Chile tsunami em signals
% Hx1  = -(Htr+j*Hti)*100*4.*pi;
% Hy1  =  (Hpr+j*Hpi)*100*4.*pi;
% Hr1  =  (Hrr+j*Hri)*100*4.*pi;

%%
Hp= zeros([361,100,396]);
prd = 1./(fftfrq(720,1/(60/(24*3600))));
prd = [4000. prd]; % The mean flow period -> 4000 days (it should have been DC but 4000 days OK)!
base_dir = 'E:\NoBackup\tsunami\Allgeyer_Modeling\resfiles\';
for i = 1 : 120,
display(['    ' sprintf('%10.5f',prd(i)) '     ' sprintf('tsu%d',i) ' tsunami25.model ' sprintf('tsu%d',i) '.source' ' 1']);
[Hx1,Hy1,Hr1] = read_chile_res([base_dir sprintf('tsu%d',i) '.res']);
Hp(i,:,:) = Hy1(:,805:1200);
end;
%Note. Compared ts synthesized from 131 harmonics to that with 200. No
% noticeable difference. hence no need of going to harmonics higher than
% 131 the saved data files are all from the first 131 harmonics
%

%%

%Hr(132:361,:,:) = 0; %filling up the higher harmonics for which the no simulations were done.
% Hp(241:301,:,:) = 0;% this is justified, since the higher harmonics of t flow are
% weaker.
Hp(121:361,:,:) = 0;

for i = 1:100,
    for j = 1:396,
%         Hr_t(:,i,j) = real(ifft([squeeze(Hr(:,i,j)).' fliplr(squeeze(Hr(2:end-1,i,j))')])); % Hr(2:end-1,i,j) because the DC to be omitted and the last value is the same (fft symm)
          Hp_t(:,i,j) = real(ifft([squeeze(Hp(:,i,j)).' fliplr(squeeze(Hp(2:end-1,i,j))')]));
%          Ht_t(:,i,j) = real(ifft([squeeze(Ht(:,i,j)).' fliplr(squeeze(Ht(2:end-1,i,j))')]));
    end;
end;

%% plot

load E:\NoBackup\tsunami\Allgeyer_Chile_simulation\time_data time_data;
fday_simulation=datenum(2010,2,27,06,34,14+time_data);
plot(fday_simulation,squeeze(Hr_t(:,50,204)));
datetick
load D:\Manoj\tsunami\2010_Chile_Easter_Island_data\ipc_data\ipc20100227vsec.sec.mat;
ipc_z = data.DATA(:,3);
fday = data.FDAY;
b1 = min(fday):(1/24):max(fday);
sp=spline(b1,ipc_z(1:100:end)'/spline(b1,eye(length(b1)),fday(1:100:end)'));
ipc_z_f=ipc_z-ppval(fday,sp);
ipc_h = data.DATA(:,1);
sp=spline(b1,ipc_h(1:100:end)'/spline(b1,eye(length(b1)),fday(1:100:end)'));
ipc_h_f=ipc_h-ppval(fday,sp);
plot(fday,ipc_z_f,'r')

%% Test
%% 
load D:\Manoj\tsunami\2010_Chile_Easter_Island_data\modeling\results_matfiles\Hr_spectra_time_series Hr_t;
load D:\Manoj\tsunami\2010_Chile_Easter_Island_data\modeling\results_matfiles\Easter_Island_ts;
%%
% Easter Islands -109.41 and -27.17
% Grid points 50,203

for i = 50:50,
    for j = 203:203,
        plot(fday,ipc_z_f,'r');
        hold on;
        plot(fday_simulation,squeeze(Hr_t(:,i,j)),'k','LineWidth' , 2);
        title(sprintf(' i = %d j = %d', i,j));
        axis([min(fday_simulation) max(fday_simulation)  -1 1]);
        pause;
        hold off;
    end;
end;


%% Movie
load E:\NoBackup\tsunami\Allgeyer_Modeling\Allegyer_Hr_ts Hr_t time_data fday_simulation;
%%
hf=figure;
kk=1;
set(hf,'Position',[680         272        1083         826]);
rect = get(hf,'Position'); 
mov = avifile('E:\NoBackup\tsunami\Allgeyer_Chile_simulation\Bz_000_4_NeetuNIO.avi','fps',5,'quality',100)
load D:\Manoj\tsunami\COLORTABLE;

for i = 100:2:500,
    h=subplot(211);
    set(h,'Position',[ 0.1 0.1 0.8 0.6]);
    AA = squeeze(Hr_t(i,:,:));
    worldmap([-39.625,-14.875],[-159,-61]);
    %geoshow(flipud(AA),[4,-14.875,-158.875],'DisplayType', 'texturemap')
    contourfm(flipud(AA),[4,-14.875,-158.875],30);
    geoshow('landareas.shp','FaceColor','none'); %
    colormap(map);
    caxis([-1 1]);
    plotm(-27.17,-109.41,'k.','MarkerSize',20,'color','green'); % Easter Islands
    textm(-27.17,-109.41+2,'IPM','FontSize',16,'color','blue'); % Easter Islands
    plotm(-35.846,-72.719,'k.','MarkerSize',20,'color','blue');%epicenter
    textm(-35.846,-72.719+2,'Epicenter','FontSize',16,'color','blue');%epicenter

    hh=colorbar('horizontal');
    set(hh,'position',[0.4007 0.1766 0.25 0.0323])
    set(hh,'FontSize',16);
    title(['B_z at sea level          ' sprintf('%d Minutes',i)],'FontSize',16);
    
    
    h1=subplot(411);
    set(h1,'position',[0.1300    0.75    0.7750    0.2]);
    plot(fday,ipc_z_f,'b');
    hold on;
    plot(fday_simulation(1:i) + 2/1440,squeeze(Hr_t(1:i,50,204)),'r','LineWidth',2);
    plot(fday_simulation(i) + 2/1440,squeeze(Hr_t(i,50,204)),'r.','MarkerSize',20);
    set(gca,'FontSize',16);
    xlabel('UTC - Hours');
    ylabel('B_z (nT)');
    legend('Observed','Predicted','location','NorthWest');
    datetick;
    axis([datenum(2010,2,27,6,0,0) datenum(2010,2,27,18,0,0) -1 1]);
    hold off;
    
     F = getframe(hf);
     mov = addframe(mov,F);
     clf;
    end;
mov = close(mov);
clear mex % clears all open avi files useful when sharing violation problem occurs  


%% plot time series of Model Vs obs

load E:\NoBackup\tsunami\Allgeyer_Modeling\Easter_Island_ts;

subplot(311)
plot(fday,ipc_z_f,'b');
set(gca,'FontSize',16);
ylabel('B_z (nT)');
hold on
plot(fday_simulation+ 5/1440,ipc_sim_z_allegyer,'r','LineWidth',2);
axis([datenum(2010,2,27,8,0,0) datenum(2010,2,27,16,0,0) -1 1]);
datetick('x','keeplimits');


subplot(312);
plot(fday,ipc_x_f,'b');
set(gca,'FontSize',16);
ylabel('B_x (nT)');
hold on
plot(fday_simulation+ 5/1440,-ipc_sim_x_allegyer,'r','LineWidth',2);
axis([datenum(2010,2,27,8,0,0) datenum(2010,2,27,16,0,0) -1 1]);
datetick('x','keeplimits');

subplot(313);
plot(fday,ipc_y_f,'b');
set(gca,'FontSize',16);
ylabel('B_y (nT)');
hold on
plot(fday_simulation+ 5/1440,-ipc_sim_y_allegyer,'r','LineWidth',2)
axis([datenum(2010,2,27,8,0,0) datenum(2010,2,27,16,0,0) -1 1]);
datetick('x','keeplimits');


%% quiver

%load coast
incr = 3;
quiver(OX(1:incr:end,1:incr:end),OY(1:incr:end,1:incr:end),U(1:incr:end,1:incr:end),V(1:incr:end,1:incr:end),5);
hold on
plot(long,lat,'k')
axis([ -140   -60   -50     -10])
set(gca,'FontSize',16);
